 Reduced.Run1<-function (y,X,TimeIndex,n=iteration,B0, burnin=burninstage,beta0,
                         ystar.start, beta.start, betafixed=TRUE, S.start,
                         P.fixed, M, likeihoodComputed=FALSE, a00, b00,
                         ystar.return=FALSE)
  {
    
    # initial values
    mcmc <- n+burnin
    N <- length(y)
    T <- length(unique(TimeIndex))
    covariate.size <- ncol(X)
    coefficient.length <- length(beta.start)*M
    if (likeihoodComputed==1){
       likelihoodave <- rep(NA, mcmc)
     }
    State <- matrix(NA, nrow=mcmc, ncol=T)
    State[1,] <- S.start
    if (ystar.return==TRUE){
      ystarm <- matrix(NA, nrow=mcmc, ncol=N)
      ystarm[1,] <- ystar.start
    }
    #loop begins
    for(g in 2:(n+burnin)) { 
      if ((g%%50 ==0)| g==n){
        cat("g=",g,"\n")
    }
      
    #######################################
    ## Blog 1: Data Augmentation
    ######################################
    if(g==2){
      state.current <- S.start
       beta.current <- matrix(beta.start, nrow=M, byrow=TRUE)
    }else{
      state.current <- State.new
      if (betafixed!=1){
         beta.current <- Beta
      }
    }

     ystar <- rep(NA, N)
     entry <- 1:N
     if (betafixed!=1){      
       Beta <- matrix(NA, ncol=covariate.size, nrow=M)
     }
     for (i in 1:M){
         group <- which(state.current==i)
         TT <- length(group)
         if (TT==0) betadrawn <- rmvnorm(1, beta0, B0)
         else{
           beta.state <- beta.current[i,]  
           XX <- data.frame(X)[TimeIndex%in%group,]
           YY <- y[TimeIndex%in%group]
           indexY <- entry[TimeIndex%in%group]
           NN <- length(YY)
           total<-0
           ystar.M <- rep(NA, NN)
           for (k in 1:NN) {
             mu <- as.matrix(XX[k,]) %*% cbind(beta.state) 
             ystar.M[k] <- BTnormlog(x=YY[k], mu=mu)
            if (betafixed!=1){
              total<-total+outerself(as.matrix(XX[k,])[1,])
            }
           }
          if (betafixed!=1){
            B1<-solve(total+solve(B0))
            beta1<-B1%*%(t(XX)%*%ystar.M+solve(B0)%*%beta0)
            betadrawn<-rmvnorm(1, beta1, B1)
          }
          ystar[indexY] <- ystar.M
         }
         if (betafixed!=1){
           Beta[i,] <- betadrawn
         }
    }

    if (ystar.return==TRUE){
      ystarm[g,] <- ystar
    }
      
    ######################################
    ## Block 2: Update the states
    ######################################
    bigP <-  trans (M, PP=P.fixed)
    filtering.out <- Filtering(Y=ystar, X=X, TimeIndex=TimeIndex, bigP=bigP,
                               Theta=beta.current, pred.return=TRUE, M=M)
    filtering.prob <- filtering.out[[1]]
    predprob <- filtering.out[[2]]
    State.new <- backward.smoothing(filter=filtering.prob, bigP=bigP,
                                    N=T, M=M)
    State[g,] <- State.new
      
    ######################################
    ## Block 3: Compute the likelihood
    ######################################
    if (likeihoodComputed==1){
       likelihood <- c()
       for (t in 1:T){
         entry <- 1:N
         group <- entry[TimeIndex==t]
         likelihoodt <- sapply(c(1:M),
         function(k){Likelihood(y=ystar[group],yo=y[group], x=X[group,],
                                theta=beta.current[k,])})
         predt <- predprob[t,]
         likelihood[t] <- sum(likelihoodt*predt)
       }
       
       likelihoodave[g] <- sum(log(likelihood))
       # cat("loglike = ", likelihoodave[g], "\n")
      }
    }
    ####################################
    ## Return MCMC output
    ####################################
    if (likeihoodComputed==1){
      return(list(likelihoodave[-c(1:burnin)], state=State[-c(1:burnin),]))
    }else{
      return(list(state=State[-c(1:burnin),], ystar=ystarm[-c(1:burnin),]))
    }
  }
